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Mechanical graphene, which is a spring-mass model with the honeycomb structure, is investigated. 
The vibration spectrum is dramatically changed by controlling only one parameter, spring tension at 
equilibrium. In the spectrum, there always exist Dirac cones at K- and K’-points. As the tension is 
modified, extra Dirac cones are created and annihilated in pairs. When the time reversal symmetry 
is broken by uniform rotation of the system, creation and annihilation of the Dirac cones result in 
a jump of the appropriately defined Chern number. Then, a flip of the propagation direction of the 
chiral edge modes takes place, which gives an experimental way to detect the topological transition. 
This is a bulk-edge correspondence of the classical system. We also demonstrate the other important 
concept, symmetry protection of the topological states, is at work in the classical system. For the 
time reversal invariant case, the topological edge modes exist for the fixed boundary condition but 
not for the free boundary condition. This contrast originates from the symmetry breaking at the 
free boundary. 

PACS numbers: 63.22.Rc,03.65.Vf,45.20.D- 


I. INTRODUCTION 

Graphene [T] , a two-dimensional crystal with a honey¬ 
comb array of carbon atoms, has been one of the most 
hot topics in condensed matter physics in this decade. 
The most prominent feature of graphene is existence of 
massless Dirac fermions, or Dirac cones at the Fermi en¬ 
ergy. In general, not limited to graphene, Dirac cones 
allow a lot of properties distinct from usual metals or 
semiconductors, such as unique responses to electronic 
and magnetic field mm, or characteristic edge states 
mm- One interesting subject of study is a manipulation 
of Dirac cones. For instance, shifting the Dirac cones in 
momentum space induces a gauge field without breaking 
the time reversal symmetry [^. By inducing a mass of 
Dirac cones (gap opening), we have a chance to observe 
topologically nontrivial phases. For instance, by appro¬ 
priately breaking the time reversal symmetry, graphene 
can go into a quantum Hall state without external mag¬ 
netic field |6]. Or, if the spin degrees of freedom and 
the spin-orbit coupling are explicitly taken account of, 
graphene becomes a time reversal invariant Z 2 topologi¬ 
cal insulator [7]. 

One direction of recent developments in Dirac and 
topological systems is exporting the concept to systems 
other than conventional solid materials. For instance, 
manipulation of Dirac cones is experimentally realized in 
an optical lattice cold atoms [8]. The other example is 
electromagnetic field in a photonic crystal, which is gov¬ 
erned by the classical Maxwell equation. As for triangu¬ 
lar and honeycomb lattices, the Dirac cones in photonic 
band are demonstrated [9HT2]. Furthermore, topologi¬ 
cally nontrivial phase is achieved in photonic crystals or 
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a coupled resonator system UMl]- Very recently, other 
classical systems, i.e., mechanical systems obeying the 
Newton’s equation of motion are also discussed in the 
context of Dirac cones or topological edge states p(6U24] . 
One great advantage of these kinds of artificial systems 
is their controllability, which enables us to access param¬ 
eter region unreachable in solids. We can easily realize 
phenomena that are difficult in solids such as merging of 
Dirac cones [75] . 

In this paper, a honeycomb spring-mass model [75|, we 
dubbed it as mechanical graphene, is investigated as a 
typical mechanical system in which Dirac physics plays a 
role. We propose a simple and feasible way, just stretch¬ 
ing the system uniformly and isotropically, to manipu¬ 
late Dirac cones. The uniform stretch modifies tension 
of springs at equilibrium, which controls the relation be¬ 
tween the transverse and longitudinal waves. Then, the 
frequency spectrum varies as a function of the tension. 
As the tension grows from zero, the frequency band struc¬ 
ture first looks like the one in p-orbital honeycomb op¬ 
tical lattice [77] (phase I), then becomes like the one in 
bilayer graphene with the trigonal warping [28| (phase 
II). Finally the spectrum is given by that of graphene 
(phase III), but with extra degeneracy. Merging of the 
Dirac cones is observed at the two critical points. This 
topological transition associated with the merging of the 
Dirac cones is detected as a jump in the appropriately 
defined Chern number [7TJ |79|, when the time reversal 
symmetry is broken by uniform rotation of the system. 
The change in the Chern number is also observed by the 
propagation direction of the “chiral edge modes” EOIISII. 
Not only proposing a practical method for Dirac cone ma¬ 
nipulation, we also demonstrate the state-of-art idea of 
topological phenomena, symmetry protection, is at work 
in this system by focusing on edge states and boundary 
condition. Without the time reversal symmetry break¬ 
ing, in-gap edge states appear for the fixed boundary 
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FIG. 1. (a) Schematic picture of the mechanical graphene, 
(b) Definitions of Ro, x', and x”. (c) The unit vectors ai 
and 02 , and the vectors connecting the nearest neighbor sites, 
(i = l,2,3). 


condition, while they do not appear for the free bound¬ 
ary condition. In the topological view point, symmetry 
breaking at the edge is behind this contrast. 


II. MODEL AND METHODS 

The model treated in this paper is composed of mass 
points, springs, and additional potential energy sources, 
where the mass points are aligned in the two-dimensional 
honeycomb lattice. [See Fig. [^a).] Motion of the mass 
points is also restricted in the two dimensional space, i.e., 
the out of plane motion is prohibited. Parameters char¬ 
acterizing this model are mass of the mass points to, a 
spring constant k, distance between the nearest neighbor 
mass points i? 0 : ^Lnd natural length of the springs Iq. In 
addition, there are some parameters to characterize ad¬ 
ditional potential energy other than the one from springs 
connecting the mass points. The mass to is fixed to I 
for simplicity from now on. rj = Rq/Iq is not necessar¬ 
ily unity, but when it deviates from unity, each string 
gives finite force even in equilibrium, which means that 
we should make the system to avoid shrinkage as a whole. 
For instance, we need to support the system by appro¬ 
priate choice of the boundary condition. Or, the po¬ 
tential source other than the springs need to be added. 
As far as the shrinkage is prevented, the mass points 
are stable at the equilibrium position since the forces 
from three springs connected to one mass point cancel 
out. Dynamical variables for this model is written as 
XRa = yRa)j where R denotes a lattice point, and 

a is a sublattice index. As we will see below, the model 
has commonality with the nearest neighbor tight-binding 
model of real graphene, but the major difference is that 
the mass point allows to oscillate in two directions, x and 
y. There is also the sublattice degrees of freedom, which 
means that the frequency spectrum for this model has 
four bands. 

For this model, elastic energy of the spring for a given 
pattern of mass point configuration has a special impor¬ 
tance. We first see it for a pair of mass points, intro¬ 
ducing parameters Rq, x' , and x” as shown in Fig. 0b)- 
Here, Rq is a vector connecting two neighboring mass 
points in equilibrium. In this paper, we adopt an as¬ 


sumption that all springs have good linearity, i.e., the 
energy of the spring Us is always written as 

Us = ^Kil-lo)^, ( 1 ) 

where I is the length of the spring at that moment, inde¬ 
pendent of the value of 1. Then, up to the second order in 
5x = x' — x”, which is required to investigate harmonic 
oscillation around the equilibrium alignment, we have 

Us = ]^i^{{Ra-lof + ‘2'{RQ-lQ)Ro-5x+5xf,-f^^^5x^) ( 2 ) 

with 

= (3) 

and Rq = iio/|.Ro|- Summation over indices appearing 
as a pair is assumed (/i, v = x,y). This formula indicates 
that when i] = 1, Us strongly depends on the angle be¬ 
tween Rq and Sx, while when rj = 0, Us is independent 
of the direction of Sx. For r; = 1, no force is given by a 
spring if Sx is normal to Rq, since such a motion only 
leads to the rotation of the spring and does not increase 
the elastic energy in the lowest order approximation. On 
the other hand, for rj — 0, there is finite force even for 
Sx = 0, and then the motion with Sx T Rq can in¬ 
crease the energy since there is a finite component of 
force normal to Rq for such motion. In short, the value 
of r] controls a ratio of strength of the restoring force for 
Sx II Rq and Sx T Rq. In other words, it controls the re¬ 
lation between the longitudinal and the transverse wave 
modes. 

Then the Lagrangian of this system becomes 

C = T-V,-V2-Vq, ( 4 ) 

with 

( 5 ) 

Ra 

^1 = ( 6 ) 
Ra 

= ^^R'a-^Rb)l%_^_RS^R'a-XRb), ( 7 ) 

{R' aRh) 

Vq = K ^%,,^_ft.S^R'a-^Rb)^ ( 8 ) 

{R' aRh) 

where Ra is a position of the sublattice a associated with 
the lattice point at R, and = {Rq — Iq)R>^/\R\. Sum¬ 
mation over {R!aRh) means that we pick up pairs of the 
nearest neighbor mass points. The term T represents 
a kinetic energy. The term Vi is contributed from the 
potential energy other than the springs connecting the 
mass points, expanded in a series of XRa up to the sec¬ 
ond order. For simplicity, we assume that this potential is 
isotropic and sublattice independent. The terms V 2 and 
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Va come from the elastic energy of the springs. As far as 
the equilibrium is achieved for xjia = 0 , the term V 3 is 
identically zero, and does not contribute to the equation 
of motion. Assuming the periodic boundary condition 
and introducing as 

^Ra = (9) 

k 

the Lagrangian is rewritten as 

^ = (10) 

k 

^ E ^ka^-ka “ ^ E KbKa^-kb^ (H) 

CL ab 

where = (r(k))af,;b,,- Here, 

fW=K + 3 m -|)) i +(. J_^j r,eW),(12) 


and 


rAB(k) = -473 + + e-‘'““^72), (13) 


with 



[See Fig. [^c) for the definitions of iig M 
Now equation to be solved becomes 


“L + Era:(fc)<6 = 0. (17) 

b 

We seek for the solution periodic in time by introducing 
a new dynamical variable ^a/x(fe) as = e“^*(/>a/i(fc)- 
Then, what we have to solve is 

-4ct>a4k) + J2Kbik)4Ak)=0, (18) 

b 

and to have a nontrivial solution, we must have 

Det[-w2i + f(fc)] =0. (19) 

Practically, this is done by diagonalizing r(fe). In that 
sense, it is possible to map the problem to a quantum 
mechanical one, by regarding r(A;) as a quantum Hamil¬ 
tonian. 

Here, we make a comment on the symmetry of r(fc). 
Except the constant diagonal elements, r(fe) has no term 


connecting a same kind of sublattice components. Since 
the constant diagonal elements only add a constant con¬ 
tribution to 4, we say that r(fc) has a “chiral” symme¬ 
try, since r(fc) anticommute with T = diag(l, 1, —1, —1) 
if the constant diagonal part is subtracted. That is, 
f'(fc) = f (fe) — (k° -I- 3k( 1 — ^))i satishes 

{f'(fe),T} = 0. (20) 

In fermionic systems, the chiral symmetry plays impor¬ 
tant roles in various situations. For instance, it stabilizes 
Dirac cones in 2D cases [5^ . 

Now, the problem is quite similar to the phonon prob¬ 
lem in graphene |33H37] . However, what is essential in the 
following arguments is to control ry in a wide range, which 
is difficult to realize in actual graphene. Furthermore, the 
effective model for graphene phonon does not respect the 
“chiral symmetry”, which plays essential role in consid¬ 
ering the topological origin of the edge modes. This is 
because if the vibrational motion in graphene is modeled 
using a spring-mass model, it requires springs connect¬ 
ing the next nearest neighbor sites and more. Lastly, our 
model neglects the motion in perpendicular to the plane 
of the honeycomb lattice. 


III. MANIPULATION OF DIRAC CONES 


A. Bulk Frequency Spectrum 


Now, let us investigate dispersion relation of uj. Before 
going to general values of ?y, we consider two limits t] = 0 
and ?7 = 1. For 7 = 0, TAsik) becomes 

TABik) = -41 + + e-‘'=“^) . (21) 


Then 4 looks exactly like a dispersion relation of the 
nearest neighbor tight-binding model for graphene (NN- 
graphene) but with extra double degeneracy. In 7 = 
0 limit, there is no distinction between the longitudinal 
and transverse wave mode, and in this case, oscillation 
in X- and j/-direction decouple, which leads to doubled 
NN-graphene dispersion. On the other hand for 7 = 1 , 
i^AB^k) becomes 



i'ABik) = — 


with 4 = k ■ ai. Then, as shown in Ref. mi the sec¬ 
ond and third bands of 4, dispersion is same as NN- 
graphene, but the first and the fourth bands are exactly 
flat and stick to the top and bottom of the other bands. 
Interestingly, r(fc) at 7 = 1 is essentially same as the 
Hamiltonian of p-orbital honeycomb optical lattice model 
introduced in Ref. EH 

Now, we follow the evolution of the dispersion relation 
as 7 changes from 1 to 0. In the following, we always 
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FIG. 2. Dispersion relations for several rj. (a) r} = 5/6. (b) rj = 2/3. (c) rj = 3/5. (d) rj 
Dirac cones in the Brillouin zone. Cyan and blue dots represent the Dirac cones. + and - 
of the inversion operator Eq. (26l. The spring constant k is scaled as k = kq/{I — ri/2). 


= 0. Insets show the positions of the 
signs on the M-point are eigenvalues 


scale the spring constant k as k = kq/(1 — ri/2) with 
kq = I, in order to make the total band width of the 
frequency dispersion constant. This scaling also makes 
the frequency of the Dirac point constant. The dispersion 
relation for rj = 5/6 is shown in Fig. [^a). It is similar 
to the one for 77 = 1 and the Dirac cone is still there at 
the K-point (and K’-point), but the first and the fourth 
bands gains finite dispersion and they are no longer flat. 
At rj = 2/3, as shown in Fig. ^b), the gap at the M-point 
in between the second and the third bands closes. In the 
vicinity of the M-point, the dispersion is linear in F-M 
direction but quadratic in M-K direction. Such a linear 
vs quadratic dispersion relation is often characteristic for 
merging of a pair of Dirac cones [53]. In fact, soon after 
77 becomes less than 2/3, there appear two new Dirac 
cones, other than those at the K- and K’- points, on M- 
K lines. [See Fig. [^c).] Due to the three fold rotational 
symmetry, there appear six new Dirac cones in the whole 
Brillouin zone in total. The new Dirac cones move from 
the M-point to the K-point as rj approaches to zero. For 
77 < 2/3, especially for small 77 where the positions of the 
new Dirac cones get close to the K- and K’-points, the 
dispersion relation looks like the one for bilayer graphene 
with the trigonal warping [55]. Finally, at 77 = 0, the 
new Dirac cones are absorbed to the originally existing 
Dirac cones at the K- and K’-points, and the dispersion 
becomes degenerated NN-graphene dispersion. 

The transition at 77 = 2/3 that is associated with the 
Dirac cone merging can be characterized by the “Herring 
number”. For electronic systems, the Herring number is 
defined as [35| 

( 23 ) 

r—1 

where d is spatial dimension, k^s are the time reversal 
invariant momenta, and N^{kr) is the number of occu¬ 
pied states with parity eigenvalue ±1 at kj.. This num¬ 
ber is also used by Fu and Kane to obtain Z 2 topological 
number of the inversion symmetric topological insulator 


with the spin-orbit coupling, or to check the existence of 
Dirac cones without the spin-orbit coupling [3^. For 2D 
cases, parity of Nh determines the parity of the num¬ 
ber of Dirac cones in half of the Brillouin zone as ex¬ 
plained below. In Ref. UHl in order to detect Dirac cones 
in generic 2D systems, we have used the Berry phase de¬ 
fined as 

m\)= E / (^nfc II I I ) 5 

where fey and k± are two momenta in independent di¬ 
rections, and is a Bloch wave function IISIIIII. 

When the system respects both of the time reversal and 
spatial inversion symmetry, 0 (fc||) is quantized into 0 or 
TT for the spinless case |3D|. As far as the quantization 
is kept, 0 (fc||) has to show a jump when it changes as a 
function of fcy , but the jump should be associated with 
a singularity of the band structure that is nothing more 
than a Dirac cone. Having this in mind, we can say that 
if 6>(0) ^ 9{'k) mod 27r, odd number of Dirac cones exist 
in the region 0 < fcy < tt, i.e., half of the Brillouin zone, 
while if 0(0) = 0 ( 7 r) mod 27r, the number of Dirac cones 
in the half Brillouin zone is even. On the other hand, the 
inversion symmetry gives us a relation |39l I42j 

4 

gi(eM-e(o)) ^ ^ (25) 

r—1 

Therefore, the Herring number has an ability to detect 
the number of Dirac cones. This idea has been applied 
to the 2D organic Dirac fermion systems |42l|43]. Now, 
we are handling a classical mechanical system, not an 
electronic system, but the Herring number is still useful. 
We have seen that the number of Dirac points in the half 
Brillouin zone is odd (one) in 77 > 2/3 and even (four) 
in 77 < 2/3. This difference should be captured by the 
Herring number. In Figs. [ll:a)-[ 2 i:d), the eigenvalues of 
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FIG. 3. Edge spectra as a function of the momentum along 
the edge for several rj — 1 with the fixed boundary condition, 
(a) rj = 5/6. (b) rj = 2/3. (c) rj = 3/5. (d) rj — 0. 



the inversion operator Uj, which is defined as 

/O 0 1 0 
V _ I 0 0 0 1 
~ 1000 
\0 1 0 0 

at the M-point are indicated for each band. We notice 
that at the transition, r] = 2/3, the positive and nega¬ 
tive eigenvalues at the M-point are interchanged. Since 
there are three distinct M-points in the Brillouin zone, 
this interchange results in the parity change of the Her¬ 
ring number. Note that Herring originally considers only 
three-dimensional cases, and in his paper, Nh is derived 
to detect the number of degeneracy loops in 3D Brillouin 
zone, instead of Dirac cones. 



B. Edge States and Symmetry Protection 


Next, we investigate edge states. In this case, we only 
apply the Fourier transformation in one direction and 
for the other direction, real space representation is used. 


Then, we obtain the equation similar to Eq. (191, but 


with larger number of components in the matrices and 
vectors. In the following, we concentrate on the zigzag 
edge. In addition, we use the fixed boundary condition, 
that is, the mass points at the boundary are connected to 
the springs that are fixed to the wall, instead of simply 
removing the springs at the boundary. [See Figs. |^b) 
and s^)-] Importance of the boundary condition is dis¬ 
cussed later. The obtained frequency spectra for several 



FIG. 4. (a) Edge spectrum for the free boundary condition 
with rj = 1. (b,c) Schematic pictnres of the (b) fixed and (c) 
free boundary conditions. 


values of 77 as functions of k \\, momentum parallel to the 
edge, are shown in Figs. [^a)-|^d). For 77 = 5/6, as in 
the case of 77 = m, we observe in-gap edge modes near 
fc|| = 0. The edge modes forms a flat band, since the 
“chiral symmetry” is preserved in this case. At 77 = 0, 
the system exactly inherits the property of NN-graphene. 
For instance, the edge modes are found near fey = tt in¬ 
stead of fc|| = 0 found for 77 > 2/3. In short, the position 
of the edge modes in the edge Brillouin zone is switched 
from near A:|| = 0 to fey = tt as 77 changes from 1 to 0. We 
observe that the new Dirac cones emerged from the M- 
points bring new edge modes near k\\ = tt and take away 
the originally existing edge modes near A:y = 0. Again, 
note that our model does not necessary leads to the re¬ 
sults exactly same as the case of the actual graphene 
phonon problems [44l |45] . 


To see the importance of the boundary condition, ui for 
77 = 1 obtained with the free boundary condition, where 
the springs at the boundary are simply removed, is shown 
in Fig. [^ Note that 77 yf 1 does not go well with the free 
boundary condition since the finite tension results in the 
deformation of the equilibrium lattice as the springs at 
the boundary becomes absent. Then, we have to take 
the deformation into account, but instead of doing that, 
we focus on 77 = 1. Figure]^ shows that there is no edge 
mode found as in-gap states. This difference between the 
two boundary conditions is deeply related to symmetry 
protection of topological phases, which is one of the most 
important concept in modern study of topological phe¬ 
nomena. In the real graphene case, it is possible to think 
a topological origin of the edge modes [1]. For that, the 
chiral symmetry plays an important role in defining the 
bulk topological number, the quantized Berry phase, and 
protecting in-gap edge states. For the fixed boundary 
condition, the “chiral symmetry”, Eq. (20), is preserved 
even after the boundary is introduced, since the springs 
connected to the mass points at the boundary is left. On 
the other hand, the free boundary condition does not re¬ 
spects the “chiral symmetry” since the absence of the 
springs at the boundary leads to the nonuniform diago- 
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FIG. 5. Real space picture of the typical eigenmodes that are localized to the edge. Upper panels: r] = 1. Lower panels: 
77 = 1/3. For rj = 1 {rj = 1/3), the neighboring mass points at the edge oscillate in-phase (anti-phase). 


nal elements in f(fc). This indicates that even the bulk 
topological property is same, the topological edge states 
are not necessarily preserved for the free boundary condi¬ 
tion, in contrast to the fixed boundary condition, which 
demonstrates the manifestation of the symmetry protec¬ 
tion. 


Now, we know that the edge modes for 77 ^ 1 and 77 ^ 0 
has different momenta along the edge. Then, the oscilla¬ 
tion pattern in the real space should look different. In or¬ 
der to convince this expectation, we perform a calculation 
using a rectangular system with long zigzag edges and 
short armchair edges. This calculation also gives a result 
to the situation that is easily accessed in experiment. In 
Fig. 1^ the typical eigenmodes whose frequency is close 
to v3, a frequency of the bulk Dirac points, is shown 
for 77 = 1/3 and 77 = 1. For 77 = 1, the mass points at 
the boundary basically oscillate in phase with the neigh¬ 
boring mass points on the edge with large wave length 
background, reflecting the fact that the edge modes are 
existing near fey = 0 in the momentum space. On the 
other hand, for 77 = 1/3, the neighboring mass points 
on the edge tend to oscillate anti-phase, reflecting the 
edge mode position in the momentum space. For 77 = 1, 
the mass points at the boundary oscillate approximately 
normal to the boundary, while for 77 = 1/3, oscillation 
in parallel to the boundary is also allowed. In Fig. 
only the typical eigenmodes are shown, but we have con¬ 
firmed that the above statements are basically valid for 
the mode near w = y/Z. Thus, the oscillation pattern 
of the eigenmode localized at the boundary can be used 
to detect the phase transition associated with the Dirac 
cone merging. This is important since the observation of 
the edge mode with finite size system may be much easier 
than measuring the dispersion relation itself in detail. 


IV. CHERN NUMBER AND CHIRAL EDGE 
MODES 

Lastly, we consider effects of rotation, which induces 
time reversal symmetry breaking. Uniform rotation as 
a whole brings two new terms in the Lagrangian, one is 
the centrifugal force and the other is the Coriolis force. 
For simplicity, the centrifugal force is neglected. The 
centrifugal force is second order in the angular frequency 
of the rotation D, while the Coriolis force is first order 
in D, and the second order term can be neglected in the 
situation that D is regarded as a small quantity. Then, 
taking account of the Coriolis force only, the equation to 
be solved becomes [2T] 

Det[-a;2i-h2ia;D-f f(fe)] =0, (27) 


where 







(28) 


In this case it is no longer possible to obtain allowed w 
by diagonalizing r(fc). Instead, it is possible to obtain 
allowed w by explicitly evaluating the determinant ana¬ 
lytically, and solving a quartic equation of uP'. When tP 
is large, the equation to be solved approximately becomes 


Dethw^i -f 2iv^D -h f(fe)] = 0, 


(29) 


and diagonalization of 2\'/hPCi + r(fc) is sufficient to ob¬ 
tain oj. Interestingly, 2\^/PPtl -|- r(fc) is exactly like the 
Hamiltonian for n-orbital honeycomb optical lattice in 

Ref. EH 

In Ref. Ell it is shown that the finite D induces a gap 
at the Dirac point for 77 = I. By solving Eq. (27), it is 


confirmed that such a gap opening remains at work as 
77 is gradually reduced from 1 . At some critical point 
77 = 77 c, the gap is closed at the M-point. Then, for 77 less 
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FIG. 6. The Chern number as a function of rj for = 0.02 
and 0.20. 


than rjc, the system is again fully gapped except rj = 0 
where the gap is again absent, rjc depends on fl, but it 
is close to 2/3 in the small 17 limit. 

After obtaining w, eigenmodes are derived as nontrivial 
solutions for 


[-uj'^i + 2iujn + t{k)]cf>nk=0 (30) 


Then, regarding the solution 4>nk as a Bloch wave func¬ 
tion \unk)i we can define the Chern number C as 


^ = E ^ / / dfcidA:: 

n=l,2 ^ 




nk 


d(j)nk dcj)^ 


nk 


\ dki dk2 


dkn 


^4^nk ^ 
dki ) 


(31) 

except at r] = 0 and rj = rjc where the gap closing occurs. 
The summation is taken over n = 1,2, since we are now 
focusing on the gap between the second and third bands. 
The Chern number can be used as a topological order 
parameter for this system. There is a well established 
method to evaluate this quantity numerically [46]. In 
Fig.je] rj dependence of the numerically obtained Chern 
number for 17 = 0.02 and 0.20 is shown. The result at 
r] = 1 is consistent with the previous work |21] . For 
small 17, the transition from Chern number 1 to —2 is 
observed near rj = 2/3. This is well understood from the 
fact that each Dirac point carries one-half contribution 
to the Chern number. For 17 = 0.2, the transition point 
deviates from 2/3, but the transition from 1 to —2 still 
exists. Transition involving change of the Chern number 
is reported in Ref. 1211 but here we propose a transition 
with completely different mechanism. 

In order to observe physical phenomena associated 
with this transition of the Chern number, we perform 
a calculation with finite system with a triangular shape 
having zigzag edges. The Chern number for the mechan¬ 
ical system itself is not a physical observable but the 
corresponding edge states reflecting nontrivial bulk can 
be observed experimentally m- Here, again, the fixed 
boundary condition is applied, i.e., the system must be 
fit in a triangular frame supporting mass points at the 
boundary by springs. We assume that the system is at 
rest, or no vibration mode is excited at t = 0. Then, 


we pick up one of the mass points and apply force that 
is sinusoidal in time with angular frequency 17/ i.e., we 
consider a forced oscillation problem. Now, equation to 
be solved can be rewritten as 

^=iX(t) + /(t) (32) 

where X{t) = *{x{t),x{t)) and 

i = (f f) ■ (33) 


Here, F is a dynamical matrix and 17 is defined as 


/no 0 
17= 6 ^0 



(34) 


f(t) denotes an external force, and is explicitly written 
as f{t) = /osinl7't. The elements of /o have a finite 
value when it corresponds to the picked mass point. A 
formal solution of this equation is 


X{t)=f dTexp(A(t - T))/(r) (35) 

Jo 


when the system is at rest for t = 0. Then, X{t) at 
any time is evaluated by diagonalizing A and perform 
the integration in Eq. (35). 

The selected snapshots of the time evolution for 77 = 
1/3 (The Chern number is —2) and 77 = 1 (The Chern 
number is 1) calculated with 17 = 0.05 are shown in 
Fig.[7| where the color map indicates the kinetic energy 
of each mass point averaged over time period 27r/17Q with 


17q = -s/S ~ 17'. In the figure, thick arrows indicates the 
mass point on which the external force is applied. Here, 
the direction of force is chosen to be normal to the edge. 
We use 17' = 1.765 for 77 = 1/3 and 17' = 1.732 for 77 = 1 
in order to excite edge modes efficiently. Owing to these 
choices, for both of 77 = 1/3 and 77 = 1, the oscillation 
amplitude is localized near the edge. As in the case of 
the rectangular system, the neighboring mass points at 
the edge oscillate anti-phase for 77 = 1/3 and in-phase for 
77 = 1. Following the time evolution of the oscillation, 
we notice that the oscillation amplitude propagates in 
the fixed direction, which indicates existence of the chi¬ 
ral edge modes. When the “wave front” of the oscillation 
reaches to the corner, it goes to the next edge instead of 
being reflected. For 77 = 1/3 with the Chern number —2, 
the wave front shows clockwise motion, while for 77 = 1 
with the Chern number 1, it shows counterclockwise mo¬ 
tion. An interesting point of this observation is that the 
difference between 77 = 1/3 and 77 = 1 resides in the 
relation between the transverse and longitudinal modes, 
and the symmetry of the system is same for 77 = 1/3 
and 77 = 1. Nevertheless, it reverses the direction of the 
motion of the wave front. 
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FIG. 7. Snapshots of the time evolution of the system. The 
color map indicates the kinetic energy of each mass point 
averaged over the time range 27r/%/3. The external force is 
applied to the mass points indicated by thick arrows. Left 
panels: t; = 1/3 and C = —2. We observe clockwise motion 
of the edge modes. Right panels: 77 = 1 and G = 1. We 
observe counterclockwise motion of the edge modes. 


Before closing, we make a comment on the choice of 
n'. In principles, the chiral edge mode is a gapless exci¬ 
tation. However, as we are handling the hnite size sys¬ 
tem, the gap is inevitably open, and we have to be careful 
on choosing H' to excite the edge modes. There are two 
kinds of finite size effects leading to gap opening. Firstly, 
the finite size effect simply makes the energy levels dis¬ 
cretized, but this does not make much problems in ex¬ 
citing the edge modes. Secondly, in some case, the edge 
modes localized on different edges are mixed by the finite 
size effect and lead to a gap that is larger than the one 
expected from the simple level discretization. In such a 
case, we have to avoid to place in such a large gap 
in order to excite the chiral edge modes effectively. For 


77 = 1/3, we encountered with the situation that such a 
gap is induced, and then, careful choice of is essential 
to observe the behavior like Fig. 

V. SUMMARY 

To summarize, it is shown that the dispersion rela¬ 
tion of the mechanical graphene is dramatically changed 
by controlling only one parameter, the spring tension at 
equilibrium. When there is no tension at equilibrium, 
the frequency dispersion is characterized by Dirac cones 
at the K- and K’-points and flat bands at the bottom 
and top of the dispersion. As the tension at equilib¬ 
rium is strengthen, extra Dirac cones are created at the 
M-points. Then, the newly created Dirac cones migrate 
from the M-points to the K- or K’-points, and the disper¬ 
sion becomes similar to the electronic dispersion of the 
bilayer graphene with the trigonal warping. Eventually, 
the extra Dirac cones are absorbed to the Dirac cones at 
the K- and K’-points, and the dispersion is characterized 
by doubly degenerated Dirac cones in the strong tension 
limit. The singular bulk dispersion is reflected in the 
edge spectrum, which is confirmed by a calculation with 
the ribbon geometry and the fixed boundary condition. 
Corresponding to the generation and merging of the bulk 
Dirac cones, the position of the edge states in the momen¬ 
tum space changes, and this change can be detected by 
observing the real space pattern of the oscillation of the 
edge modes. On the other hand, the free boundary con¬ 
dition gives no edge modes. In the topological viewpoint, 
absence of the edge mode for the free boundary condition 
is explained by breaking of the “chiral symmetry” at the 
edge. These observation forms a good example of the 
bulk-edge correspondence. Lastly, Coriolis force is intro¬ 
duced to discuss the time reversal symmetry breaking. 
Then, generation and merging of the Dirac cones found 
without Coriolis force appear as a jump in the Chern 
number, which is defined regarding the vector of nor¬ 
mal modes as Bloch wave functions. Furthermore, it is 
demonstrated that when the sign of the Chern number is 
flipped, then the propagation direction of the chiral edge 
modes is really flipped. Importantly, flip in the propa¬ 
gation direction is caused by only controlling the tension 
at equilibrium, which keeps the symmetry of the system 
intact. That is, the system looks same, but still the edge 
modes propagate in the opposite directions. 
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